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ABSTRACT 

This  paper  gives  guidelines  for  the  development  of  computer  programs  for 
the  numerical  simulation  of  semiconductor  devices.  For  this  purpose  the  basic 
mathematical  results  on  the  corresponding  elliptic  boundary  value  problem  are 
reviewed.  In  particular,  existence,  smoothness  and  structure  of  the  solutions 
of  the  fundamental  semiconductor  equations  are  discussed.  Various  feasible 
approaches  to  the  numerical  solution  of  the  semiconductor  equations  are  described. 
Much  emphasis  is  placed  on  constructive  remarks  to  help  authors  of  device  simula¬ 
tion  programs  to  make  decisions  on  their  code  design  problems.  In  particular, 
criteria  for  an  optimal  mesh  generation  strategy  are  given.  The  iterative  solution 
of  the  systems  of  nonlinear  and  linear  equations  obtained  by  discretising  the  semi¬ 
conductor  equations  is  discussed.  An  example  is  given  showing  the  power  of  these 

concepts  combined  with  modem  numerical  methods  in  comparison  to  classical  approaches 
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SIGNIFICANCE  AND  EXPLANATION 


Many  different  codes  for  the  simulation  of  semiconductor  devices  such  as 
transitors,  diodes,  thyristors  are  already  circulated.  Most  of  them  solve  the 
basic  set  of  semiconductor  equations  in  the  steady  state  case,  which  represents 
a  nonlinear  system  of  three  second  order  elliptic  equations.  During  the  last 
15  years  this  problem  has  also  been  the  subject  of  analytical  investigations  by 
researchers  from  different  disciplines.  Ibis  paper  reviews  how  some  results  of 
these  investigations  can  be  used  to  improve  the  performance  of  numerical  methods 
for  solving  the  semiconductor  equations.  Hie  qualitative  analysis  of  the  problem 
shows  how  appropriate  finite  difference  and  finite  element  methods  can  be  con¬ 
structed  and  what  criteria  have  to  be  used  in  an  adaptive  mesh  selection  strategy 
in  order  to  require  a  minimal  amount  of  gridpoints  while  still  providing  a 
sufficiently  accurate  solution.  Various  questions  concerning  the  solution  of  the 
large,  sparse,  nonlinear  system  of  algebraic  equations  which  arise  in  these  calcu 
lations  are  also  discussed. 
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1.  Introduction 

The  characteristic  feature  of  early  device  modeling  is  the 
separation  of  the  interior  of  the  device  into  different  regions, 
the  treatment  of  which  could  be  simplified  by  various  assumptions 
like  special  doping  profiles,  complete  depletion  and 
quasineutrality.  These  separately  treated  regions  were  simply 
put  together  to  produce  the  overall  solution.  If  results  in  an 
analytically  closed  form  are  intended,  any  other  approach  is 
prohibitive.  Fully  numerical  modeling  based  on  partial 
differential  equations  /61/  which  describe  all  different  regions 
of  semiconductor  devices  in  one  unified  manner  was  first 
suggested  by  Gummel  /29/  for  the  one  dimensional  bipolar 
transistor.  This  approach  was  further  developed  and  applied  to 
pn-junction  theory  by  De  Mari  /13/,  /14/  and  to  IMPATT  diodes  by 
Scharfetter  and  Gummel  /SO/. 

A  two  dimensional  numerical  analysis  of  a  semiconductor 
device  was  carried  out  first  by  Kennedy  and  O'Brien  /35/  who 
investigated  the  junction  field  effect  transistor.  Since  then 
two  dimensional  modeling  has  been  applied  to  fairly  all  important 
semiconductor  devices.  There  are  so  many  papers  of  excellent 
repute  that  it  would  be  unfair  to  cite  only  a  few.  Recently  also 
the  first  results  on  three  dimensional  device  modeling  have  been 
published.  Time  dependence  has  been  investigated  by  e.g.  /37/, 

/ 44/  and  models  in  three  space  dimensions  have  been  announced  by 
e.g.  /8/ ,  /ll/,  /67/ ,  /68/ . 
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Zn  spite  of  all  these  important  and  successful  activities, 
the  need  for  economic  and  highly  user  oriented  computer  programs 
became  more  and  more  apparent  in  the  field  of  device  modeling. 
Especially  for  MOS  devices  which  have  evolved  since  their 
invention  by  Kahng  and  Atalla  /32/  to  an  incredible  standard, 
modeling  in  two  space  dimensions  has  become  inherently  important 
because  current  flow  controlled  by  a  perpendicular  field  is  an 
intrinsically  two  dimensional  problem.  One  such  program  which 
has  been  applied  successfully  in  many  laboratories  is  called 
CADDET  /59/.  We  have  also  tried  to  bridge  that  gap  and  developed 
MIN I MOS  /53/,  /51/  for  the  two  dimensional  static  analysis  of 
planar  MOS  transistors. 

2.  Analvsis  of  the  Static  Semiconductor  Equations 


In  this  chapter  we  review  some  of  the  existing  analytical 
results  for  the  fundamental  semiconductor  equations  concerning 
existence  and  structure  of  their  solutions.  These  results  are  of 
importance  in  both  the  theoretical  and  practical  context,  since  - 
as  we  will  see  in  the  next  chapter  -  the  knowledge  of  the 
structure  and  smoothness  properties  of  solutions  is  indeed 
essential  for  the  development  of  a  numerical  solution  method. 
The  most  familiar  model  of  carrier  transport  in  a  semiconductor 
device  has  been  proposed  by  Van  Roosbroeck  /61/.  It  consists  of 
Poisson's  equation  (2.1),  the  current  continuity  equations  for 
electrons  (2.2)  and  holes  (2.3)  and  the  current  relations  for 


These  relations  form  a  system  of  coupled  partial 
differential  equations.  Poisson's  equation,  coming  from 
Maxwell's  laws,  describes  the  charge  distribution  in  the  interior 
of  a  semiconductor  device.  The  balance  of  sinks  and  sources  for 
electron-  and  hole  currents  is  characterized  by  the  continuity 
equations.  The  current  relations  describe  the  absolute  value, 
direction  and  orientation  of  electron-  and  hole  currents.  The 
continuity  equations  and  the  current  relations  can  be  derived 
from  Boltzmann's  equation  by  not  at  all  trivial  means.  It  is  not 
our  intention  to  present  in  this  paper  the  ideas  behind  these 
considerations.  The  interested  reader  is  refered  to  /61/  and  its 
secondary  literature  or  text  books  on  semiconductor  physics  e.g. 
/7/,  /31/ ,  /52/ ,  /5«/. 


2.1  The  Validity  of  the  Basic  Semiconductor  Equations 

Xt  is  of  prime  importance  to  be  aware  that  equations  (2.4) 
and  (2.S)  are  not  capable  to  describe  exactly  all  phenomena 
occur ing  in  real  devices.  For  instance,  they  do  not  characterize 
effects  which  are  caused  by  degenerate  semiconductors  (e.g.  heavy 
doping) .  /38/,  /60/,  /63/  discuss  some  modifications  of  the 
current  relations,  which  partially  take  into  account  the 
consequences  introduced  by  degenerate  semiconductors  (e.g. 
invalidity  of  Boltzmann's  statistics,  bandgap  narrowing).  These 
modifications  are  not  at  all  simple  and  lead  to  problems 
especially  in  the  formulation  of  boundary  conditions  /47/,  /62/. 
Xn  case  of  modeling  .  MOS  devices,  degeneracy,  owing  to  the 
relatively  low  doping  in  the  channel  region,  is  practically 
irrelevant.  For  modern  bipolar  devices,  though,  bearing  in  mind 
shallow  and  extraordinarily  heavily  doped  emitters,  it  is  an 
absolute  necessity  to  account  for  local  degeneracy  of  the 
semiconductor. 


Just  as  further  examples  (2.4)  and  (2.5)  do  not  describe 
velocity  overshoot  phenomena  which  become  apparent  at  feature 
lengths  of  O.ljhn  for  silicon  and  ljfcn  for  gallium-arsenide  /25/. 


Certainly  no  effects  which  are  due  to  ballistic  transport  (the 
existence  of  which  is  still  questionable  /30/)  are  included.  The 
latter  start  to  become  important  for  feature  sizes  below  O.OlJhn 
for  silicon  and  O.ljhn  for  gallium-arsenide  /26/.  Considering  the 
state  of  the  art  of  device  miniaturization,  neither  effect  has  to 
bother  the  modelists  of  silicon  devices.  For  gallium-arsenide 
devices  new  ideas  are  mandatory  in  the  near  future  /25/ ,  /46/, 
/  45/ . 


2.2  Domain  and  Boundary  Conditions 

Host  of  the  existing  programs  which  solve  the  semiconductor 
equations  are  restricted  to  a  rectangular  device  geometry.  This 
is  not  essential  as  far  as  the  analysis  of  the  equations  is 
concerned.  In  this  chapter  we  shall  assume  that  the  equations 
(2.1)-(2.5)  are  posed  in  a  domain  D  of  Rn  (n=l,2,3)  with  a 
piecewise  smooth  boundary  3d.  Equations  (2.1)-(2.5)  are  subject 
to  a  mixed  set  of  Oirichlet  and  Neumann  boundary  conditions. 
That  means  3d  consists  of  three  parts  3d=3DjU3d2i/3d3  .  dDj^ 
denotes  the  part  of  the  boundary  where  the  device  is  surrounded 
by  insulating  material.  There  one  assumes  the  boundary 
conditions: 

3^/3 nj_«  3n/3nJ_*  9p/3nJ_»  0  (2.6) 

Here  nj_  denotes  the  unit  normal  vector  on  3d  which  exists 
anywhere  except  at  a  finite  number  of  points  (arbitrarily  defined 
corners  of  the  simulation  geometry) .  3d2  denotes  the  part  of  the 

boundary  corresponding  to  the  ohmic  contacts.  There  V,  n  and  p 
are  prescribed.  The  boundary  conditions  can  be  derived  from  the 
applied  bias  %  and  the  assumptions  of  thermal  equilibrium  and 
vanishing  space  charge: 

*•%  +  %uilt-in'  ■  "i’-  -  -  P  -  =  -  0 

The  last  two  conditions  in  (2.7)  can  be  rewritten  as: 


^  ,  ▼  .■»  .-v  ”  L  *  /•  i.  1  .  *  _  •  r  * ' 


n  -  (^  C2+4-ni2  +  C)/2 

p  -  (^|  C2+4-ni2  -  C)/2 


(2.8) 


In  many  applications  it  is  desired  to  consider  controlled 
insulator-semiconductor  interfaces  (e.g.  MOS  devices).  So  3d3 
denotes  the  part  of  the  boundary  which  corresponds  to  such  an 
interface.  There  we  have  the  interface  conditions: 


S„"«l  "  ^p'n|  “  0 

1 sem  'ins 


(2.9) 


Again  nl  denotes  the  normal  vector  on  3d.  C _  and  €■  _ 

*  sem  ins 

denote  the  permittivity  constants  for  the  semiconductor  and  the 

insulator  respectively.  3v/3n||  and  3M»/3nh  denote  the 

•sem  'ins 

onesided  limits  of  the  derivatives  perpendicular  to  the  interface 
approaching  the  interface.  Within  the  insulator  the  Laplace 
equation:  div  grad  q» •  0  holds. 

2.3  Dependent  Variables 

For  analytical  purposes  it  is  often  useful  to  use  other 
variables  than  n  and  p  to  describe  the  system  (2.1)-(2.5).  Two 
other  sets  of  variables  which  are  frequently  employed  are 
(V*Vn»Vp)  and  (M*,u,v)  which  relate  to  the  set  (*p,n,p)  by: 


-  n^e^^t,  P  -  ni.e(^/Ut 

-  n^e^t-u,  p  -  ni-e'W°t-v 


(2.10) 

(2.11) 


(2.10)  can  be  physically  interpreted  as  the  application  of 
Boltzmann  statistics.  However  (2.10)  a. so  can  be  regarded  as  a 
purely  mathematical  change  of  variables  so  that  the  question  of 


the  validity  of  the  Boltzmann  statistics  does  not  need  to  be 
considered.  The  use  of  (V,^,^)  a  priori  excludes  negative 
carrier  densities  n  and  p,  which  may  be  present  as  undesired 
nonphysical  solutions  of  (2.1)-{2.5)  if  we  use  (tp,n,p)  or  (tp,u,v) 
as  dependent  variables.  As  we  will  see  later  in  this  chapter  the 
advantage  of  the  set  (q*,u,v)  is  that  the  continuity  equations 
(2.2)/  (2.3)  and  current  relations  (2.4),  (2.5)  become 
self-adjoint.  This  also  has  an  important  impact  on  the  use  of 
iterative  schemes  for  the  solution  of  the  evolving  linear  systems 
(cf.  chapter  4).  However,  owing  to  the  enormous  range  of  the 
values  of  u  and  v,  the  sets  (M>,n,p)  or  (H»,Vn»Vp)  have  to  be 
prefered  for  actual  computations.  We  personally  favour  the  set 
(S»,n,p)  . 


2.4  The  Existence  of  Solutions  and  Scaling 


The  basic  answer  to  the  question  of  existence  of  solutions 
can  be  found  in  Mock  /43/  or  under  slightly  different  assumptions 
in  Bank,  Jerome  and  Rose  /S/.  Both  proofs  are  based  on 
Schauder's  fixpoint  theorem.  They  are  both  valid  for  arbitrarily 
shaped  domains  and  boundary  conditions  of  the  type  previously 
described  without  an  interface  (9d^*{}).  Both  papers  consider 
the  case  of  vanishing  generation/recombination  rate  (R-0  in 
(2.2),  (2.3)).  In  the  setting  of  Mock  (t|>,u,v)  is  used  as 
dependent  variables.  The  equations  are  scaled  so  that  the 
intrinsic  carrier  density  n^,  the  thermal  voltage  and  the 
ratio  elementary  charge/permittivity  are  equal  to  unity.  Thus, 
combining  the  continuity  equations  (2.2),  (2.3)  and  current 
relations  (2.4),  (2.5),  we  have  the  system: 


div  grad 


w 

eT*u 


-V 

e  *  v  - 


(2.12) 


div  (e^grad  u)  *  0 


(2.13) 


div 


grad  v) 


0 


(2.14) 
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Then  a  map  M:tp+  y  ia  defined  (details  in  /43/  or  /4/)  such 
that  the  evaluation  of  M  requires  the  solution  of  (2.13)  and 
(2.14)  and  a  fixpoint  of  M  (M(<p*)-MM  together  with  the 
according  functions  (u.v)  is  a  solution  of  the  whole  system 
(2.12)-(2.14) .  The  existence  of  a  fixpoint  is  shown  by 
Schauder's  fixpoint  theorem.  Questions  concerning  the  degree  of 
smoothness  of  these  solutions  (the  existence  of  derivatives)  are 
discussed  in  /42/. 

However.  Schauder's  theorem  is  not  constructive  and  does  not 
indicate  that  iterating  the  map  M  will  actually  lead  to  the 
fixpoint.  Moreover,  it  does  not  give  any  information  about  the 
structure  of  the  solution  which  is  of  vital  interest  for  actual 
computations.  Since  the  dependent  variables  in  the  system 
(2.1)-(2.5)  are  of  different  order  of  magnitude  and  show  a 
strongly  different  behaviour  in  regions  with  small  and  large 
space  charge  the  first  step  towards  a  structural  analysis  of 
(2.1)-(2.5)  has  to  be  an  appropriate  scaling.  A  standard  way  of 
scaling  (2.1)-(2.5)  has  been  given  by  De  Mari  /14/.  There  is 
scaled  by  the  thermal  voltage  Ut,  n  and  p  are  scaled  by  n^ 
(similar  to  Mock  /43/)  and  the  independent  variables  are  scaled 
such  that  all  multipying  constants  in  Poisson's  equation  become 
unity.  Although  physically  reasonable  this  approach  has  the 
disadvantage  that  n  and  p  in  general  are  still  several  orders  of 
magnitude  larger  than  <|».  A  scaling  which  reduces  tp,  n  and  p  to 
the  same  order  of  magnitude  has  been  given  by  Vasiliev'a  and 
Butuzov  /65/.  This  approach  makes  the  system  (2. l)-(2. 5) 
accessible  to  an  asymptotic  analysis  which  is  given  together  with 
applications  in  /40/.  /41/  and  /39/.  There  n  and  p  are  scaled  by 
the  maximum  absolute  value  of  the  net  doping  C  and  the 
independent  variables  are  scaled  by  the  characteristic  length  of 
the  device.  More  precisely  the  following  scaling  factors  are 
employed. 


quantity 

symbol 

value 

.  •*  *.  ^  a  , 

X 

1 

max(x-y) ,  x.y  in  D 

* 

°t 

k*T/q 

(2.15) 

n.p 

4 

max | C | 
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After  scaling  the  equations  become: 


X2*div  grad  4»  *  n  -  p  -  C  (2.16) 

diw  (  grad  n  -  n*grad  )  *  -R 
div  {  grad  p  +  p*grad  4»  )  “  "R 

Here,  for  simplicity  only,  J>n  and  have  been  assumed  to  be 
constant.  It  should  be  noted  that  the  following  analysis  also 
holds  if  the  usual  smooth  dependence  of  and  P  on  n,  p  and 
grad  e.g./54/  is  assumed.  Since  the  independent  variable  x  has 
been  scaled,  equations  (2.16)  are  now  posed  on  a  domain  D  with 
maximal  diameter  equal  to  one.  The  small  constant  X2  multiplying 
the  Laplacian  in  (2,16)  is  the  minimal  Debye  length  of  the 
device: 


X2  -  ------  (2.17) 

l2-q.«f 

1  and  «t  are  defined  in  (2.15).  Thus  for  high  doping  «f>>l) 
X2  will  be  small.  For  instance  for  a  silicon  device  with 
characteristic  length  25Jhn  and  «l*1020cm-3  we  compute  for  X2  at 
approximate  room  temperature  T=300K:  X2*4.10~^°. 


R  denotes  again  the  scaled  generation/recombination  rate. 
In  the  analysis  given  in  /41/  the  usual  Shockley-Read-Hall  term 
has  been  used  which  after  scaling  is  of  the  form: 


R  « 


n  +  p  +  2-  nrX)2 


r-  1/2 


(2.18) 


R  is  in  general  a  (not  necessarily  mildly)  nonlinear 
function  of  n,p  and  grad*p.  Thus  different  models  of  R  may 
influence  the  analytical  results  quite  drastically.  This  is 
obviously  to  be  expected  as  in  many  operating  conditions  the 
device  behaviour  depends  strongly  on  the  net 
generation/recombination  R. 


2.5  The  Singular  Perturbation  Approach 


(2.16)  represents  a  singularly  perturbed  elliptic  system 
with  perturbation  parameter  X.  The  advantage  of  this 
interpretation  is  that  we  can  now  obtain  information  about  the 
structure  of  solutions  of  (2.16)  by  using  asymptotic  expansions: 
In  the  subdomains  of  Ds  where  the  solutions  behave  smoothly  we 
expand  them  into  power  series  of  the  form: 

oo 

w(x,A)  •Zv~i(x)‘Xi,  w=(4»,n,p)T  (2.19) 
i-0 

which  implies  a  smooth  dependence  on  X.  C  -  the  scaled  doping  - 
is  smooth  in  these  subdomains  and  exhibits  a  sharp  transition 
across  tr.e  pn-junctions  in  the  device.  For  the  case  of  an  abrupt 
junction  this  behaviour  is  represented  by  a  discontinuity  across 
an  n-1  dimensional  manifold  P:(x=x(s),  s  of  Rn_1)  in  the  device. 
Thus  r  is  a  point  in  1  dimension,  a  curve  in  2  dimensions  and  a 
surface  in  3  dimensions.  Of  course  one  curve  or  surface  has  to 
be  used  for  each  junction.  Since  the  procedure  is  the  same  for 
each  of  the  junctions  it  is  demonstrated  only  for  one  junction. 
In  the  case  of  an  exponentially  graded  doping  profile  C  consists 
of  two  parts: 

C  -  C"  +  C*  (2.20) 

where  C”and  C*  are  discontinous,  C"  is  piecewise  smooth  and  C~  is 
exponentially  decaying  to  zero  away  from  P.  In  the  vicinity  of  P 
the  expansion  (2.19)  is  not  valid  and  has  to  be  supplemented  by  a 
"layer”  term  acording  to  the  singular  perturbation  analysis: 


w(x,A) 


■  E[wJ(x)  +  wj  (s,  t/X) )  •X1', 
i*0 


w»(q>,n,p)T 


(2.21) 
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Here  the  following  coordinate  transformation  has  been 
employed:  For  a  point  in  the  vicinity  of  P  s  denotes  the 
parameter  value  at  the  nearest  point  on  P  and  t  denotes  its 
distance  perpendicular  to  P  (cf.  Fig.  1).  Thus  the  solution  of 
the  semiconductor  equations  exhibits  internal  layers  at 
pn- junctions. 

The  wj  and  w£  in  (2.21)  can  now  be  determined  separately  and 
the  structure  of  the  solution  is  given  by  its  partition  into  the 
smooth  part  ZwT-X1  and  its  rapidly  varying  part  Ewt'X1.  w"  has 
to  satisfy  the  reduced  equations: 


0 


(2.22) 


div 

(grad  n~  -  n~-grad<lT)  =  -R“ 

(2.23) 

div 

(grad  p"  +  p"*grad4T)  =  -R~ 

(2.24) 

For  the  sake  of  simplicity  but  without  loss  of  generality 
the  mobilities  p^  and  P have  been  assumed  to  be  constant. 
(2.22)— (2.24)  is  subject  to  the  boundary  conditions  (2.6)-(2.9) . 
Of  course  the  condition  of  vanishing  space  charge  is  redundant 
with  (2.22).  Since  C"  is  discontinous  at  P  and  (2 . 22) - (2 . 24) 
represents  a  second  order  system  of  two  equations  four  "interface 


conditions"  have 

to  be  imposed  at  P. 

They  are  of  the  form: 

o  1 x=x- 

.  •  Ve-%|.  * 

1  x=x+ 

(2.25) 

eo  ' x=x- 

“  Po*e,^U  ^ 

'x-x+ 

(2.26) 

a. 

3  J  ^  t  1 

n°’ nll x=x+ 

(2.27) 

^o-Sllj.5. 

*  J  1  1 

Po*nll^5?+ 

(2.28) 

and  w|^+  denote  the  onesided  limits  of  w  as  x  tends  to 


where  w|-_ 

P  from  each  side,  nl  denotes  the  unit  normal  vector  on  P.  J“ 

*  n 


and  3“  are  the  zeroth  order  terms  of  the 
Po 

(scaled)  electron  and  hole  current  densities 


smooth  parts  of  the 
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(2.29) 


-  grad  no  “  nc*grad  qr 
J"o  -  grad  pS  +  pj-grad  qr 

(2.22) -(2.24)  together  with  (2.25)-(2.28)  and  the  boundary 
conditions  (2.6)-(2.9)  define  the  reduced  problem  whose  solution 
is  an  0(A)  approximation  to  the  full  solution  away  from  f1.  As  we 
will  see  in  the  next  chapter  the  reduced  problem  is  a  useful  tool 
for  the  development  and  analysis  of  numerical  methods,  since  it 
(especially  the  conditions  (2.25)-(2.28) )  has  to  be  solved 
implicitly  by  any  discretisation  method  which  requires  a 
reasonable  number  of  grid  points. 

The  equations  for  the  rapidly  varying  parts  w£  reduce  to 
ordinary  differential  equations.  That  means  that  only 
derivatives  with  respect  to  the  "fast*  variable  t/A  occur.  Since 
the  rate  of  decay  of  wj  depends  heavily  on  q»  the  width  of  the 
layer  grows  with  the  applied  voltage)  a  fact  which  is  absolutely 
well  known  by  device  physicists,  but  which  becomes  nicely 
apparent  by  the  singular  perturbation  approach. 

3.  Numerical  Solution  of  the  Semiconductor  Equations 

In  this  chapter  we  discuss  some  of  the  problems  occur ing  in 
the  numerical  solution  of  the  semiconductor  equations  and  the 
analysis  of  existing  numerical  methods.  From  the  viewpoint  of 
numerical  analysis  there  are  essentially  four  major  topics  to  be 
considered.  The  first  one  is  the  type  of  discretisation  to  be 
used.  There  exist  programs  for  both  Finite  Element  and  Finite 
Difference  discretisations  of  the  system  (2.1)-(2.5).  As 
outlined  in  the  previous  chapter  the  solution  exhibits  a  smooth 
behaviour  in  some  subregions  of  the  domain  whereas  in  others  it 
varies  rapidly.  Thus  a  nonuniform  mesh  is  mandatory  and  adaptive 
mesh  refinement  is  desirable.  So  the  second  topic  is  the 
question  how  to  set  up  the  mesh  refinement  algorithm  i.e.  which 
quantities  have  to  be  used  to  control  the  mesh.  Each  type  of 


discretisation  will  lead  to  a  large  sparse  system  of  nonlinear 
equations  and  so  the  solution  of  this  system  is  the  third  topic. 
As  fourth  topic  we  discuss  linear  equations  solvers  which  have  to 
be  used  in  topic  three.  For  topics  one  to  three  many  methods 
have  been  designed  especially  for  the  semiconductor  equations. 
These  points  will  be  discussed  in  this  chapter.  For  topic  four 
standard  numerical  analysis  is  commonly  used  and  so  its 
discussion  will  be  deferred  to  chapter  four.  For  the  sake  of 
simplicity  in  nomenclature  we  shall  only  consider  the 
two-dimensional  case  in  this  chapter.  However,  all  results  given 
in  the  following  can  be  generalized  to  three  dimensions  in  a 
straightforward  manner.  So,  the  equations  are  posed  in  a  domain 
D  of  |R  and  x  *  (x,y)  denotes  the  independent  variable. 

3.1  Discretisation  Schemes 

Using  Finite  Elements  or  Finite  Differences  one  has  to  take 
into  account  that  Poisson's  equation  (2.1)  is  of  a  different  type 
than  the  continuity  equations.  Poisson's  equation  -  in  the 
scaling  of  Markowich  /40/  using  the  variables  (4»,u,v) 

X2-div  grad  q>  *  e^-u  -  e"^v  -  C  (3.1) 

is  a  singularly  perturbed  elliptic  problem  whose  right  hand  side 
has  a  positive  derivative  with  respect  to  Thus  it  is  of  a 
standard  form  (as  discussed  in  e.g.  /22/)  except  for  the 
discontinous  or  exponentially  graded  term  C.  Equations  of  that 
type  are  generally  well  behaved  and  it  suffices  to  apply  a  usual 
discretisation  scheme.  In  the  case  of  Finite  Differences 
equation  (3.1)  is  discretized  by: 

X2 •  (div  gradhm  ^  -  n^  -  p^  -  Cfx^y^)  (3.2) 
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Ei+l/2,j  *  <*i+l,j-*i,j>/hi 


(3.3) 


Ey 

Ei,j+l/2 


ri,j+i 


hi  *  xi+i-xi'  Kj  " 

(dlv  grad  Wifi  -  2-(Ej+1/2  j  -  eJ_1/2, ^/Oy^i-l* 

♦  J-<ET.jn/2  *  *(,i-i/2,/<kjrti-i’ 


(3.4) 


Here  q^,  nij  and  Pij  denote  the  approximations  to  q»,  n  and 
p  at  the  gridpoint  (x^.yj).  E*+i/2,j  denotes  the  value  of  3q*/3x 

at  (xi+i/2“<xi+xi+l)^2'  yjJ*  Ei,j+l/2  denotes  the  value  of  3q»/3y 
at  (x^,  y j+i/2* (y j+1) /2) .  If  one  of  the  neighbouring 
gridpoints  (xi+1,yj)#  (xi_1,yj),  (x^y^),  <xi'yj-i>  does  not 
exist  -  as  possible  in  a  terminating  line  approach  /!/,  /2/  or  in 
the  Finite  Boxes  approach  /24/  -  (3.4)  has  to  be  modified.  He 
will  go  into  some  detail  concerning  these  modifications  in  the 
next  section.  In  the  case  of  Finite  Elements  classical  shape 
functions  can  be  used  (i.e.  linear  shape  functions  for  triangular 
elements,  bilinear  shape  functions  for  rectangular  elements) . 

It  turns  out  that  the  discretisation  of  the  continuity 
equations  is  more  crucial  than  the  discretisation  of  Poissons’ s 
equation.  The  usual  error  analysis  of  discretisation  methods 
provides  an  error  estimate  of  the  form: 


max  |wh-w|  <■  c«H 


(3.5) 


wh  denotes  the  numerical  approximation  to  w(x,y)  *  (4»,n,p). 
H  denotes  the  maximal  gridspacing.  The  constant  c  will  in 
general  depend  on  the  higher  order  derivatives  of  w.  The 
singular  perturbation  analysis  /41/  shows  that  derivatives  of  qf , 
n~  and  p~  in  (2.21)  are  of  magnitude  0(A”3)  -  0(A-4)  locally  near 
the  junction  (A  is  defined  in  (2.17)).  /41/  shows  also  that, 
even  if  a  nonuniform  mesh  is  used,  the  amount  of  gridpoints 
required  to  equidistr ibute  the  error  term  in  (3.5)  can  be 
proportional  to  A  2  which  is  of  course  prohibitive.  Therefore  a 


discretisation  scheme  is  needed  where  the  constant  c  in  (3.5) 
does  not  depend  on  the  higher  derivatives  of  the  rapidly  varying 
terms  qF,  n~  and  p~.  For  the  case  of  Finite  Differences  such  a 
scheme  was  given  by  Schar fetter  and  Gummel  /50/.  They 
approximate: 


Jn  ■  grad  n  -  n*grad  q> 

div  J_  •  3j*/9x  +  JjJ'/iy  ■  R 
n  n  n 


(3.6) 


(3.7) 


JSi+1/2,j  “^<<,*'i+1#j-M»i,j)/2)*(nl+lfj-ni>j)/hi  - 
-  (ni ^  j+ni+1  #  j )  /2 *  <Vi+1 1  f  j)  /*i 

(3.8) 

Jf>i,j+l/2  "  (Vi,  j}/2) '  (ni,  j+l“ni,  j)/kj 

“  <ni,j+ni,j+l)/2' (Tl,j+l”Vi,j>/kj 
f(s)  ■  s*coth (s) 

2'  (J*i+l/2,  j  ‘  J*i-l/2,  j1/<hi+hi-x>  + 

+  2* <J"i» j+1/2  "  J"i, j-l/2,/(kj+kj-l>  *  Ri, j  (3,9) 

Jni+i/2,j  denotes  the  value  of  J*  at  (xi+1/2“ (*i+xi+i> /2 ' 
Yj)*  J«if j+1/2  denote3  the  value  of  J&  at  (x*, 
yj+l/2“(yj+yj+l)/2) *  Th*  continuity  equation  for  holes  is 
discretized  analogously.  Scharfetter  and  Gummel  gave  a  physical 
reasoning  for  the  derivation  of  their  scheme.  Markowich  et  al. 
/41/  proved  that  in  one  dimension  the  Scharfetter-Gummel  scheme 
is  uniformly  convergent.  That  means  that  the  error  constant  c  in 
(3.5)  does  not  depend  on  the  derivatives  of  qf ,  n'  and  p~  in 
(2.21)  and  therefore  not  on  X.  For  two  dimensions  /41/  shows 
that  the  choice  f(s)  -  s*coth(s)  is  necessary  for  uniform 
convergence.  Exponentially  fitted  schemes  like  the  Scharfetter- 
Gummel  scheme  have  been  analyzed  by  Kellog  /34/f  /33/  and  Doolan 
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/17/  ( for  different  classes  of  problems) .  The  reason  for  the 
uniform  convergence  of  these  schemes  is  that  inside  the 
pn- junction  layers  the  interface  conditions  (2.25)  and  (2.26)  are 
satisfied  automatically  if  |gradt|i|  is  large  and  the  gridspacing 
is  not  0 (X) . 

The  results  for  Finite  Difference  schemes  suggest  that  a 
similiar  approach  (like  the  exponentially  fitted  schemes)  should 
be  used  in  the  case  of  Finite  Elements.  This  fact  has  been 
intuitively  observed  by  Engel  /21/  for  the  one-dimensional  case. 
A  modeling  group  at  IBM  has  tried  to  make  use  of  the  Scharfetter- 
Gummel  scheme  fdr  Finite  Elements  in  two  and  three  space 
dimensions  / 9 /,  /8/,  /12/.  However,  we  have  the  impression  that 
their  approach  needs  still  quite  a  bit  of  analysis,  although  it 
has  been  used  effectively  by  other  modelists  too  e.g.  /49/. 
Macheck  /36/  has  tried  to  develop  a  more  rigorous  discretisation 
for  Finite  Elements  using  exponentially  fitted  shape  functions. 
He  uses  classical  bilinear  shape  functions  for  ¥  and 

«l(x,y)  -  [1  -  *i(x,y)l-[l  -  *2<**y>l  (3-11) 

«*2<*»y>  ■  *l(x.y)  •U-*2(x,y)] 

<3(x.y)  ■  ^(x,y)  •  *a(x,y) 

*4<*/y>  ■  U  -  *i(x,y)l*  *2(x,y) 

for  u,  and 

Pl(x,y)  ■  II  -  *1 (x,y) ) • II  -  *2(x,y)]  (3.12) 

02<*'y>  *  »l(x,y)  *11  -  #2(**y>i 

03(x,y)  -  ?l(x,y)  •  »2<}‘*y) 

^4(x,y)  ■  U  -  *j.(x,y)]*  92  (x,y) 


for  v,  where 


(3.13) 


(Pi (x,y)  ■  f(x,|^ 

(P2(x,y)  *  f(y,|^ 

»l(x,y)  -  f(x,-  |j) 

#2(x»y)  *  *(y.-  f^) 

with:  f(x,a)  »  (exp (ax) -1)/ (exp (a) -1) 


(3.14) 


The  advantage  of  these  shape  functions  is  that  they 

accomodate  nicely  the  layer  behaviour  of  the  solution.  They 
degenerate  into  the  ordinary  bilinear  shape  functions  when  the 
electric  potential  is  constant.  In  order  to  be  able  to  switch 
from  coarse  to  fine  grid  spacing  in  different  subdomains 
transition  elements  have  to  be  used  (as  outlined  in  the  next 
section).  However,  no  theoretical  investigations  have  been 

carried  out  so  far  to  analyse  the  uniform  convergence  properties 
of  this  method. 

3.2  Grid  Construction 


Since  subregions  of  strong  variation  of  n  and  p  alternate 
with  regions  where  these  quantities  behave  smoothly  (i.e.  their 
gradients  are  small)  different  meshsizes  are  mandatory  in  these 
subregions.  Thus  the  discretisation  scheme  should  be  able  to 
switch  locally  from  a  coarser  to  a  finer  grid.  For  the 
exponentially  fitted  (Scharfetter-Gummel)  Finite  Difference 
discretisation  schemes  this  is  done  by  the  Finite  Boxes  approach 
/24/.  Grid  lines  can  terminate  when  the  mesh  is  likely  to  be 
coarsened  (cf.  Fig. 2).  The  point  (x^+^»yj)  does  not  belong  to 
the  mesh.  Thus  the  equations  for  the  point  (x^,y^)  have  to  be 

modified  since  t|*+1  j#  ni+],  j  and  ?i+i  j  are  not  available.  This 
is  done  by  proper  interpolation  between  the  (j-l)-st  and  (j+l)-st 
y-level.  So  (div  grad  q*)^  is  approximated  by: 
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(div  grad  V  j  ■ 

•  2*<<kj-l*Ei+l/2,j+l  +  kj'Ei+l/2,j-l)/(kj+kj-l) 

-  ^-l^j^Wl*  + 

+  2*(Ei,j+l/2  ”  Ei,j-l/2)/(kj+kj-l) 


(3.15) 


i-1/2, j'  i, j+1/2  '‘•v"  “ 

continuity  equations  are  approximated  by 


etc.  are  defined  in  (3.3).  The 


2*<<kj-l,Jni+i/2,j+i  +  kj*Jni+l/2,j-l)/(kj+kj-l) 
“  ^1-1/2, j,/(hi*hi-l)  + 

♦  2*<Jni#j+i/2  "  Jni,j-l/2)/(kj+kj-l>  "  Ri,j 


(3.16) 


JJJi_l/2  j»  Jni  j+l/2  ete*  ar*  defined  (3.8).  For  reasons 
of  numerical  stability  only  one  gridline  is  allowed  to  terminate 
at  a  box.  This  approach  is  a  generalisation  of  the  "Terminating 
Line"  approach  introduced  by  Adler  /!/ ,  /2/  as  already  mentioned. 

Zn  the  Finite  Element  approach  of  Macheck  /36/  transition 
elements  composed  of  three  triangles  are  used  to  coarsen  the  mesh 
locally  (cf.  Fig. 3).  Within  these  triangles  a  different  set  of 
shape  functions  has  to  be  used.  They  are  derived  by  holding  the 
current  densities  3n  and  constant  along  the  edges  of  a 
triangle  similar  to  the  approach  of  /10/. 

Zn  the  Finite  Element  as  well  as  in  the  Finite  Difference 
(Boxes)  approach  the  question  arises  which  criteria  should  be 
used  to  generate  the  mesh.  If  the  user  of  a  simulation  program 
has  to  define  his  elements  or  nodes  a  priori  as  input  parameters, 
this  could  perhaps  be  done  by  experience  /10/.  However,  if  -  as 
it  is  the  case  for  modern  user  oriented  programs  -  an  adaptive 
mesh  selection  is  desired  mathematically  formulated  criteria  are 
a  "sine  qua  non".  Generally  such  criteria  should  satisfy  two 
conditions.  Firstly  they  should  not  cause  the  program  to 
construct  more  gr idpoints/elements  than  necessary  to  achieve  a 
certain  accuracy.  Secondly  they  should  guarantee  that  a 


prescribed  relative  accuracy  4  is  really  achieved  once  they  are 
satisfied.  A  usual  way  to  design  adaptive  mesh  refinement 
procedures  is  to  equidistr ibute  the  local  truncation  error  of  the 
discretisation  scheme.  In  the  case  of  Finite  Differences  this 
error  is  proportional  to  the  meshsize  and  the  third  and  fourth 
derivatives  of  tp,  n  and  p.  Markovich  /41/  however  showed  that  it 
is  practically  not  possible  to  equidistr ibute  this  quantity.  In 
the  case  of  a  simple  MOS-transistor  0 (4~2X~2)  gridpoints  would  be 
required.  On  the  other  hand  the  singular  perturbation  analysis 
shows  that  the  solution  of  the  difference  scheme  approximates  the 
solution  of  the  reduced  problem  (2. 22) - (2, 24)  even  if  this 
criterion  is  not  satisfied  inside  the  layer  regions  (inversion 
layer  and  space  charge  regions) .  Therefore  the  quantity  to  be 
equidistr ibuted  is  the  discretisation  error  of  Poisson's  equation 
(i.e.  the  partial  derivatives  of  the  space  charge  times  the 
meshsizes) .  This  equidistr ibution  can  be  relaxed  inside  the 
pn-junctlon  layers  by  e.g.  simply  limiting  the  number  of 
gridpoints  there. 

3.3  Linearisation  Schemes 

Each  discretisation  scheme  (Finite  Differences  or  Finite 
Elements)  will  lead  to  a  large  sparse  system  of  nonlinear 
equations  to  be  solved.  The  theory  of  iterative  methods  to  solve 
these  equations  is  to  a  large  extent  independent  of  the  used 
discretisation  and  so  it  is  convenient  to  view  the  whole  problem 
as  solving  a  nonlinear  system  of  equations  iteratively  by  solving 
linear  systems.  The  existing  numerical  methods  can  essentially 
be  divided  into  two  classes:  The  first  approach,  a  block 
nonlinear  iteration  algorithm,  is  due  to  Gummel  /29/  and  uses  the 
fact  that  the  current  relations  are  linear  in  the  variables  u  and 
v  (as  defined  in  (2.11)).  In  these  variables  the  equations 
become  (again  we  use  the  scaling  of  /36/) : 


(3.17) 


A2* div  grad  “  e^*u  -  -  C 

dlv  Jft  »  R,  Jn  -  e^grad  u  (3.18) 

div  Jp  ■  -R,  Jp  ■  -e”^*grad  v  (3.19) 

Gummels  approach  works  as  follows:  Given  (l|»,u,v)k  is 

computed  by  solving: 

A2*div  grad  »|^+1  -  1*uk  -  e”^+1  •vk  -  C  (3.20) 

k+1 

subject  to  the  appropriate  boundary  conditions.  Then  u 
k+1 

and  v  are  computed  from: 


div  Jk+1  «  R(grad  q^+1,uk,vk) ,  Jk+1  -  e^+1*grad  uk+1 

k4.r  (3.22) 

div  Jk+1  -  -R(grad  qft+1,uk,vk)  ,  Jk+l  -  -e"1*^  *grad  vk+1 


together  with  the  boundary  conditions  for  u  and  v.  (3.21) 

k+1  k+1 

and  (3.22)  are  two  decoupled  linear  equations  for  u  and  v  . 
Poissons' s  equation  (3.20)  is  nonlinear  in  this  setting  and 
therefore  it  has  to  be  solved  iteratively  itself  in  each  step  by 
a  Newton  like  method.  Since  Newton's  method  is  an  inner 
iteration  within  the  overall  iteration  process  (3.20)-(3.22)  it 
may  not  be  necessary  to  let  this  inner  iteration  "fully  converge” 
/27/.  It  could  for  instance  be  considered  to  do  only  one  Newton 
step  for  each  iteration.  This  would  lead  to  the  linear  equation: 

A2, div  grad  t^+1  -  (e^*uk  +  e"^-vk)  •  (^+1-<^)  + 

+  e^-uk  -  e~^-vk  -  C  (3.23) 

instead  of  (3.20).  The  advantage  of  Gummels's  method  is  obvious. 
(3.20)-(3.22)  can  be  solved  sequentially  which  decreases  the 
required  amount  of  storage  and  computing  time  drastically  for 
each  step.  However,  bad  convergence  properties  can  be  observed 
in  the  case  of  high  currents.  This  is  explained  by  viewing 


-19- 


(3. 20) - (3 . 22 )  as  iterating  the  map  M:  (uk,vk)+(uk+1,vk+1)  where 

the  evaluation  of  M  involves  the  solution  of  (3.20).  Then  the 

norm  of  the  linearisation  of  M  (as  an  operator  acting  in  the 

*  *  *  * 

appropriate  spaces)  at  the  fixpoint  M(u  ,v  ) « ( u  ,v  )  is 
proportional  to  the  current  densities  /42/. 

The  second  approach  to  the  solution  of  the  nonlinear 
equations  (2.1) -(2.5)  is  a  damped  modified  Newton  method.  To 

I* 

solve  the  general  equation  F(x)«0  one  computes  the  sequence  <x  > 
by: 

Mk*rfk  -  -F (xk) ,  xk+1  «  xk  +  tk-rfk  (3.24) 

For  the  usual  Newton  method  Mk  ■  F' (xk)  and  tk  ■  1  holds. 

Bank  and  Rose  /4/  have  given  criteria  for  the  choice  of  the 

If 

damping  parameters  t  which  guarantee  global  convergence. 

•  If 

Moreover  they  investigate  how  well  4  has  to  approximate  the 
classical  Newton  step  in  order  to  get  a  certain  rate  of 
convergence.  They  obtain  that  the  rate  of  convergence  is  p 
(l<p<2)  if: 

|Mk-rfk  +  F (xk) |  -  0(|F(xk)|P)  (3.25) 

holds  asymptotically  for  k  ♦  oo.  Alternatively  Bank  and  Rose  /3/ 
suggested  Mk  *  AkI  +  F' (xk)  where  Ak  is  proportional  to  |F(xk)|. 
Franz  /24/  tested  this  method  with  good  success.  However,  he 

I# 

additionally  chooses  damping  parameters  t  according  to  Oeuflhard 
/15/,  /16/ . 

Since  this  approach  has  the  disadvantage  that  all  three 
equations  are  solved  simultaneously  -  and  therefore  the  storage 
requirements  are  fairly  large  -  we  suggest  a  Block-Newton-SOR 
method  /24/.  Defining  F»(Fi,F2»F3)T  Newton's  method  at  step  k 
is: 
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(3.26) 
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Under  the  assumption  that  the  Jacobian  is  definite  one  can 
use  a  classical  block  iteration  scheme  (iteration  index  m)  for 
the  solution  of  the  k-th  Newton  step: 

(3.27) 
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Since  the  coefficient  matrix  of  (3.27)  is  block  lower 
triangular  one  can  decouple  the  elimination  process  into  three 
linear  systems  (3.28) - (3. 30)  which  have  to  be  solved 
sequentially. 

(3.28) 

•  ^cmfl  ■  *F]>(Sfc,nk,pk)  - 

(3.29) 

.  Ajkm+1  ■  -f2  (H^n11,^)  -  •  <^aiH-l  -  -j-2.  ^p*®1 

(3.30) 

flp,k  .  3p,k  4r>,k 

■  “F3(Mfcfnkfpk)  -  jjjj-  •  ^Jo^l  -  j —  • 
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This  iteration  method  has  (like  Gummel’s  method)  the 
advantage  that  the  equations  can  be  solved  sequentially.  To  end 
up  with  the  Block-Newton-SOR  method  one  has  to  cesubstitute  the 
series  expansions  on  the  right  hand  side  of  (3.28)-(3.30)  and  to 
introduce  a  relaxation  parameter  Ufc 


(3.31) 


.  ^|totl  m  ~qi  .  f,  jkm^pk^jJati) 

•V 


(3.32) 


j-2  .  ^nkm+l  .  -u,  .  F2(q*^^1,nk,pltWptan) 


3F,k 


(3.33) 


3  .  jyaiH-1  M  -ui  •  F3(«|K^|)ffl,+i,n,tWn,a,H'^,p*<) 


This  method  converges  linearly  /48/.  However,  we  still  have 
to  perform  thorough  investigations  in  order  to  properly  judge  the 
convergence  properties. 


4.  Solution  of  Linear  Systems 

For  any  of  the  linearization  procedures  which  have  been 
outlined  in  the  last  chapter  a  large  sparse  linear  equation 
system  (4.1)  has  to  be  solved  repeatedly. 

A*x  *  b  (4.1) 

A  has  been  derived  by  linearizing  discretized  PDEs.  Hence  A 
has  only  five  to  nine  nonzero  entries  per  row  and  block  (the 
blocks  are  defined  in  (3.26))?  A  is  very  sparse.  For  the 
solution  of  these  special  types. of  linear  systems  of  equations 
two  classes  of  methods,  can,  in  principle,  be  used:  direct 
methods  which  are  based  on  elimination  and  iterative  methods.  An 


excellent  survey  on  that  subject  has  been  published  recently  by 
Duff  /18/.  Classical  Gaussian  elimination  is  not  feasibly  for 
our  systems  of  equations  because  the  rank  of  A  in  (4.1)  is  very 
large  and  A  has  many  coefficients  which  are  zero.  Therefore, 
modifications  of  the  classical  Gaussian  elimination  algorithm 
have  to  be  introduced  to  account  for  the  zero  entries.  There 
exist  quite  a  few  activities  on  that  subject  (c.f.  /19/)  and 
powerful  algorithms  which  treat  the  nonzero  coefficients  only  are 
available  (so  called  sparse  matrix  codes) .  Another  serious 
drawback  of  direct  methods  lies  in  the  fact  that  the  upper 
triangular  matrix  which  is  created  by  the  elimination  process  has 
to  be  stored  for  back  substitution.  This  matrix  has  usually  more 
nonzero  entries  than  the  matrix  A.  Therefore,  memory  requirement 
of  direct  methods  is  substantial.  One  advantage  of  the  linear 
systems  obtained  from  the  discretised  semiconductor  equations  is 
that  no  pivoting  in  order  to  maintain  numerical  stability  is 
needed.  In  spite  of  all  drawbacks  of  direct  methods,  their  major 
advantage  is  high  accuracy  of  the  solution.  However,  we  feel 
that  for  the  semiconductor  problems  iterative  algorithms  are  to 
emphasize.  Nevertheless  we  and  many  others  have  observed 
difficulties  with  respect  to  the  convergence  speed  of  iterative 
methods,  so  that  the  direct  methods,  which  require  an  exactly 
predictable  amount  of  computer  resources,  will  always  stay  in 
consideration. 

The  fundamental  idea  of  relaxation  methods  (which  are  the 
best  established  iterative  methods)  is  the  splitting  of  the 
coefficient  matrix  A  (4.1)  into  three  matrices  D,  E,  F  (4.2). 

A  «  D  -  E  -  F  (4.2) 

D  denotes  the  diagonal  entries  of  A;  -E  denotes  a  lower 
triangular  matrix  which  consists  of  all  sub-diagonal  entries  of 
A;  and  -F  denotes  an  upper  triangular  matrix  which  consists  of 
all  super-diagonal  entries  of  A. 

With  an  arbitrary  non  singular  matrix  B  which  has  the  same 
rank  as  A  the  linear  system  (4.1)  can  be  rewritten  to  (4.3): 
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B*x  +  (A-B)  -x  >=  b  (4.3) 

One  obtains  an  iterative  scheme  by  setting: 

B*xk+1  «  b  -  (A-B) -xk  (4.4) 

k+1 

(4,4)  can  be  solved  for  x  : 

xk+1  *  (I-B_1-A)-xk  +  B-1*b  (4.5) 

The  scheme  (4.5)  will  converge  if  condition  (4.6)  holds: 

f(I-B_1-A)  <  1  (4.6) 


(4.6)  is  a  necessary  and  sufficient  condition  where  $ 

denotes  the  spectral  radius  /64/.  Any  relaxation  method  can  be 

derived  by  differently  choosing  the  matrix  B  from  the  splitting 
of  A  (4.2).  The  simplest  scheme,  the  point-Jacobi  method,  uses  D 
for  B.  Matrix  D  is  a  diagonal  matrix  and,  therefore,  is  easily 
invertible.  The  Gauss-Seidel  method  uses  D-E  for  B.  The  matrix 
D-E  is  a  lower  triangular  matrix.  Therefore  one  has  only  to 

perform  a  forward  substitution  process  for  its  inversion.  The 

successive  overrelaxation  method  (SOR)  uses  a  parameter  U>  within 
the  range  ]0,2[.  The  iteration  matrix  B  is  defined: 

B  »  O/ttl  -  E  •  (4.7) 

Since  B  is  again  a  lower  triangular  matrix,  its  inversion  is 
instantly  reduced  to  a  substitution. 

The  major  advantage  of  these  iterative  methods  lies  in  their 
simplicity.  They  are  very  easy  to  program  and  demand  only  low 
memory  requirement.  As  already  noted,  they  converge  if  condition 
(4.6)  holds.  However,  this  is  generally  difficult  to  prove.  A 
sufficient  condition  for  convergence  is  that  A  is  positive 
definite  (4.8)  which  is  the  normal  case  for  five-point-star 
discretized  PDEs. 

xT-A*x  >  0  for  all  x^O  (4.8) 
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It  should  be  noted  again  here  that  the  current  relations  and 
continuity  equations  are  not  self  adjoint  if  (4»,n,p)  are  used  as 
variables  (see  (2.10),  (2.11)).  However,  the  transformation: 


tl) 

n  •  eT*u,  p 


e  T*v 


(4.9) 


results  in  a  similarity  transformation  of  the  iteration  matrix  in 
(4.6).  Thus  the  spectral  radius  of  the  iteration  matrix  is  not 
influenced  and  the  same  convergence  properties  are  obtained  as  if 
the  system  had  been  discretized  in  its  self  adjoint  form  with 
(tp,u,v)  as  variables. 


Some  point-iterative  schemes  can  by  accelerated  quite 
remarkably  with  the  conjugate  gradient  method  or  the  Chebyshev 
method.  An  excellent  survey  on  these  topics  can  be  found  in 
/28/. 


Various  activities  can  be  observed  for  the  development  of 
more  powerful  algorithms  with  the  advantages  of  iterative 
schemes.  One  of  the  best  known  algorithms  which  has  been 
established  in  semiconductor  device  analysis  is  Stone's  strongly 
implicit  procedure  /58/.  Stone's  idea  was  to  modify  the  original 
coefficient  matrix  A  by  adding  a  matrix  N  (whose  norm  is  much 
smaller  than  the  norm  of  A)  so  that  a  factorization  of  (A+N) 
involves  less  computational  effort  than  the  standard 
decomposition  of  A.  Assuming  this  has  been  done,  the  development 
of  an  iterative  procedure  is  then  fairly  straightforward  because 
the  equation  can  be  written  as: 

(A+N) *x  «  (A+N) *x  +  (b-A*x)  (4.10) 
which  suggests  the  iterative  procedure: 

(A+N)*xk+1  »  (A+N)-xk  +  (b-A*xk)  (4.11) 

When  the  right  hand  side  is  known  and  if  (A+N)  can  be 
factorized  easily,  (4.11)  gives  an  efficient  method  for  directly 
solving  for  x^+^.  Furthermore,  one  would  intuitively  expect  a 
rapid  rate  of  convergence  if  N  is  sufficiently  small  compared  to 
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A.  We  will  refrain  from  explaining  in  detail  Stone's  suggestion 
of  how  to  choose  the  perturbation  matrix  N  because  this  has  been 
done  thoroughly  in  many  publications  e.g.  /23/,  /55/,  /58/.  A 
major  disadvantage  of  Stone's  method  is  that  it  is  only 
applicable  for  linear  systems  obtained  by  a  classical  Finite 
Difference  discretisation.  It  is  not  applicable  for  systems 
obtained  by  the  Finite  Boxes  approach  or  the  general  Finite 
Element  approach. 

There  exist  a  few  algorithms  which  are  similar  to  Stone's 
method  in  terms  of  underlying  ideas.  The  most  attractive  are  the 
method  of  Dupont  et  al.  /20/,  the  "alternating  direction 
implicit"  methods  e.g.  /6/,  /23/,  /66/  and  the  Fourier  methods 
/57/,  /64/.  However,  most  of  these  sophisticated  algorithms  lack 
general  applicability. 

No  matter  which  iterative  method  is  used  one  has  to  deal 
with  the  question  of  an  appropriate  termination  (convergence) 
criterion.  Usually  (4.12)  is  applied  with  a  properly  chosen 
relative  accuracy  €: 

|xk+1-xk|  <  C-|xk+1|  (4.12) 

Since  increments  still  accumulate  when  (4.12)  is  already 
satisfied  we  suggest  to  use  (4.13)  instead  of  (4.12): 

|xk+1-xk|  <  C-|xk+1|- (1-?{G))  (4.13) 

$  (G)  can  be  estimated  as  |xk+1-xk|/|xk-xk-1|. 

One  disadvantage  of  all  strongly  implicit  methods  and  also 
the  direct  methods  is  that  they  cannot  be  implemented  efficiently 
on  a  computer  with  a  pipe-line  architecture  (vector  processor) . 
Some  comments  on  that  subject  have  been  given  in  /18/. 


-26- 


5.  A  Glimpse  on  Results 

As  an  illustrative  example  a  relatively  simple  structure,  a 
two  dimensional  diode,  is  chosen.  Fig.  4  shows  the  doping 

1 A  —3 

profile  as  birds-eye-view  plot.  A  substrate  with  10  cm 
acceptor  concentration  and  an  exponentially  graded  n-region  with 

1  Q  .1 

10  cm  maximum  doping  is  assumed.  The  initial  mesh  is 
automatically  generated  from  the  doping  profile  and  the  geometry 
definition.  The  simulation  domain  (device  geometry)  is  a  square 
of  100JI  times  lOOJfcn  size.  At  the  n-region  an  ohmic  contact  with 
length  20J»m  is  assumed.  The  substrate  is  fully  contacted.  The 
initial  mesh  for  a  Finite  Boxes  program  is  shown  in  Fig.  5  and 
for  a  Finite  Element  program  in  Fig.  6.  The  point  allocation  is 
identical  for  both  representations.  The  grid  consists  of  121 
points  versus  178  when  all  gridlines  are  extended  throughout  the 
device.  This  clearly  demonstrates  the  advantage  of  the  Finite 
Boxes  approach.  In  Finite  Element  representation  one  has  to  deal 
with  80  rectangular  elements  and  17  transition  elements  which 
consist  of  51  triangles. 

Fig.  7  shows  the  final  grid  for  an  operating  condition  of 
0.7V  forward  bias  in  Finite  Boxes  representation.  This  mesh  is 
obtained  after  several  adaption  processes  using  the  criteria 
given  in  chapter  3.  It  consists  of  270  points  (versus  480  for 
the  classical  approach).  In  Fig.  8  the  potential  distribution  is 
drawn.  From  this  plot  and  even  better  from  the  electron  density 
(Fig.  9)  one  nicely  can  deduce  the  effects  of  high  injection. 
E.g.  the  substrate  is  flooded  with  carriers.  Fig.  10  shows  the 
magnitude  of  the  electron  current  density.  The  peak  value  is 
about  180  A/cm  .  The  sharply  pronounced  peak  which  exists  at  the 
transition  of  the  Dirichlet  boundary  condition  to  the  Neumann 
boundary  condition  corresponds  to  a  singularity  of  the  carrier 
densities.  Physically  interpreted  this  effect  is  well  known  as 
con  tac  t-cor  ne  r-current-c  r owd i ng . 


Fig.  -11  shows  the  final  grid  for  an  operating  condition  of 
-20V  (reverse)  bias  in  Finite  Element  representation.  This  mesh 


consists  of  363  points  (625  for  classical  Finite  Differences) 
which  correspond  to  277  rectangular  elements  and  41  transition 
elements  (123  triangles) .  The  electron  density  for  this 
operating  point  is  given  in  Fig.  12.  One  nicely  observes  the 
depletion  region  and  the  typical  shape  of  the  drop  of  the 
electron  density  in  that  region  owing  to  thermal  generation.  In 
Fig.  13  the  magnitude  of  the  electron  current  density  is  drawn. 
The  singularity  at  the  contact  corner  is,  although  it  still 
exists,  not  so  pronounced.  Note  that  there  are  about  seven 
orders  of  magnitude  difference  in  the  peak  value  compared  to  Fig. 
10. 

6.  Conclusion 

In  this  paper  we  have  presented  an  analysis  of  the  steady 
state  semiconductor  equations  and  the  impact  of  this  analysis  on 
the  design  of  device  simulation  programs.  By  appropriate  scaling 
we  have  transformed  the  semiconductor  equations  into  a  singularly 
perturbed  elliptic  system  with  nonsmooth  data.  Information 
obtained  from  the  singular  perturbation  analysis  has  been  used  to 
investigate  stability  and  convergence  of  discretisation  schemes 
with  particular  emphasis  on  the  adaptive  construction  of 
efficient  grids.  We  have  reviewed  algorithms  for  the  solution  of 
nonlinear  and  linear  systems  of  the  discretized  semiconductor 
equations.  An  example  has  demonstrated  the  power  and  flexibility 
a  device  simulation  program  can  achieve  when  using  the 
information  we  have  presented  for  program  design. 
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Fig.  9  Electron  concentration  (0.7V)  [cm-3]  (log) 

Fig.  10  Electron  Current  Density  (0.7V)  [A/cm2]  (lin) 

Fig.  11  Final  Mesh  for  20V  Reverse  Bias  (Finite  Elements) 
Fig.  12  Electron  Concentration  (-20V)  [cm-3]  (log) 

Fig.  13  Electron  Current  Density  (-20V)  [A/cm2]  (lin) 
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